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Abstract 

We performed 3D numerical simulations of the merger of equal-mass binary neu- 
tron stars in full general relativity using a new large scale supercomputer. We 
take the typical grid size as (505,505,253) for {x,y,z) and the maximum grid 
size as (633,633,317). These grid numbers enable us to put the outer bound- 
aries of the computational domain near the local wave zone and hence to calcu- 
late gravitational waveforms of good accuracy (within ~ 10% error) for the first 
time. To model neutron stars, we adopt a I^-law equation of state in the form 
P — [r — l)pe, where P, p, e and F are the pressure, rest mass density, specific 
internal energy and adiabatic constant. It is found that gravitational waves in 
the merger stage have characteristic features that reflect the formed objects. In 
the case that a massive, transient neutron star is formed, its quasi-periodic oscil- 
lations are excited for a long duration, and this property is reflected clearly by 
the quasi-periodic nature of waveforms and the energy luminosity. In the case 
of black hole formation, the waveform and energy luminosity are likely damped 
after a short merger stage. However, a quasi-periodic oscillation can still be seen 
for a certain duration, because an oscillating transient massive object is formed 
during the merger. This duration depends strongly on the initial compactness of 
neutron stars and is reflected in the Fourier spectrum of gravitational waves. To 
confirm our results and to calibrate the accuracy of gravitational waveforms, we 
carried out a wide variety of test simulations, changing the resolution and size of 
the computational domain. 



§1. Introduction 



Binary neutron stars such as the Hulse- Taylor binary pulsar adiabatically in- 
spiral as a result of the radiation reaction of gravitational waves, and they eventually 
merge. The latest statistical study suggests that such mergers could occur approx- 
imately once per year within a distance of about 30 Mpc in the most optimistic 
scenario. "^^ Even in the most conservative scenario, the event rate would be ap- 
proximately once per year within a distance of about 400 Mpc. This implies that 
the merger of binary neutron stars is one promising source for kilo-meter-size laser 
interferometric detectors, such as LIGO, TAMA, GEO600 and VIRGO. 

The waveform in the merger stage*^ of binary neutron stars is believed to have 

*' Hereafter, we refer to the stage after which the hydrodynamic interaction between two neutron 
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two characteristic frequencies. One is that of quasi-normal modes (QNMs) of a black 
hole that is likely formed in the final stage of the merger in most cases. Perturbation 
studies have revealed that the frequency of the fundamental quadrupole QNM, /qnM; 
is between 0.1(7rMBH)~^ and 0.5(7rMBH) ""^^ where Mbh is the gravitational 
mass of the black hole. The frequency depends strongly on the angular momentum 
parameter of the black hole q = J/Afg^ where J denotes the angular momentum, 
and is higher for larger q. (For example, for q = 0, /qnm = 0.1(7rMBH)~^ and 
for g ^ 1, /qnm ^ 0.5(7rMBH)~"'^-) Studies of quasiequilibrium configurations just 
before merger suggest that the value of q for a black hole formed after a merger is 
~ 0.8 — 1^) (see also Table I). Then, the frequency of dominant QNMs is^) 

/qnm ~ (0.25 - 0.5)(7rMBH)~' ~ (5 - 10) kHz|^^^, (1-1) 

where Mq denotes the solar mass. This frequency will be rather high, so that it may 
be difficult to detect QNMs of black holes, even with advanced laser interferometric 
detectors such as LIGO II or resonant mass detectors, in the near future. 

The other characteristic frequency is associated with fundamental oscillation 
modes of a merged massive object formed transiently after the onset of a merger. 
Numerical simulations in the frameworks of Newtonian,^-*' post-Newtonian, " -""^^ 
semi-relativistic, and fully general relativistic (GR) gravity have indicated 

that the frequency is approximately between ~ 2 and ~ 3 kHz, depending on the 
equations of state and the compactness of neutron stars. Therefore, it is also too high 
for first generation laser interferometric detectors to detect. However, in contrast 
to the frequency of QNMs of a formed black hole, it is not extremely high, and it 
may be in the frequency range for advanced laser interferometric detectors, such as 
LIGO II, or resonant mass detectors in the near future. This implies that from 
gravitational waves induced by such fundamental oscillations, it may be possible 
to obtain information on the dynamical merger stage and the nature of transient 
massive neutron stars, which will be useful in understanding the nature of neutron 
stars. 

Numerical hydro dynamic simulations employing full general relativity provide 
the best approach for understanding the waveforms from the merger of binary neu- 
tron stars. Over the last few years, numerical methods for solving coupled Einstein 
equations and hydrodynamic equations have been developed, and now such 

simulations are feasible. The purpose of this paper is to investigate the waveforms in 
the merger stage using a fully GR simulation carried out on a large scale supercom- 
puter which has become available recently. This is the first paper that investigates 
gravitational waveforms from the merger of binary neutron stars using full general 
relativity with good accuracy. 

A minimum grid number for extracting gravitational waveforms in numerical 
relativity is estimated in the following manner. Since it is impossible to carry out a 
simulation from the late inspiraling stage to the merger stage, due to the limitation 

stars sets in as the 'merger stage'. The stage before the merger stage is referred to as the 'early 
stage'. 
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on computational resources, we usually start simulations from innermost orbits of 
binaries. Assuming that a binary is in a circular orbit with a Kepler angular velocity, 
the wavelength of gravitational waves in this orbit is estimated as 

where i7io, M, and a denote the angular velocity of the innermost orbit, gravitational 
mass and orbital separation, respectively. Since the system becomes more compact 
after the onset of the merger, the typical wavelength of gravitational waves is shorter 
in the merger stage. Thus, the wavelength given in Eq. (1-2) is the longest one 
throughout a simulation. This implies that the location of the outer boundaries in 
the computational domain, L, has to be larger than Aio, i.e., 

L > Aio. (1-3) 

On the other hand, the grid size Ax has to be small enough to resolve neutron stars 
and a black hole that is likely formed after the merger in many cases. The shortest 
length scale in the simulation appears to be the radius of the formed black hole, 2M. 
Here we assume that the mass of the black hole is approximately equal to the initial 
total gravitational mass M. Requiring the radius to be covered by at least 10 grid 
points, we stipulate a upper bound as 

Ax < 0.2M. (1-4) 

From the constraints (1-3) and (1-4), the required minimum grid number for covering 
one dimension is estimated as 

where the factor 2 appears because both plus and minus directions need to be covered 
as the computational domain. 

In previous works, we performed simulations using FACOM VPP300 in the 

data processing center of the National Astronomical Observatory of Japan (NAO J) . 
The available memory of this machine was ~ 30 GBytes, so that the maximum grid 
size we used was (293,293,147) for {x,y,z). With this grid size, L was typically 
~ Ao/3, where Aq denotes the wavelength of gravitational waves for the I = \m\ = 2 
mode at t = 0. Thus, we were not able to calculate gravitational waves in the 
early stage accurately nor to investigate the convergence of waveforms while varying 
L and Ax for a wide range in the previous simulations. Since the wavelength of 
gravitational waves becomes short during the merger, L becomes comparable to the 
wavelength A, and thus the situation is improved in late stages of the merger. Even 
so, the outer boundaries were still located in the local wave zone (i.e., L ^ X), and 
therefore the accuracy of extracted gravitational waves is suspect. 

Since April 2001, a more powerful vector-parallel supercomputer FACOM VPP5000 
has been available at NAO J. About 700 GBytes memory is available on this machine. 
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Thus, it is possible to use 2-3 times as large A'' as that in previous simulations car- 
ried out on VPP300. This allows us to take L as ^ Aq, keeping a grid spacing of 
Ax ~ 0.2M. This implies that gravitational waves can be computed fairly accurately 
throughout the simulation. In addition, with this new supercomputer, it is possible 
to investigate the convergence of waveforms as the size of the computational domain 
L (and the location of the wave extraction at ~ L) is increased from ^ Ao to ~ Aq. 
In the merger stage, the wavelength of gravitational waves becomes between ~ Ao/3 
and ~ Ao/2, since a merged object is compact. As a result, the accuracy of the ex- 
tracted gravitational waveforms is significantly improved, and it is feasible to draw 
scientifically meaningful conclusions regarding gravitational waveforms, luminosity 
and spectra. 

The paper is organized as follows. In §2, wc review the basic equations, gauge 
conditions and methods for setting initial conditions that we currently adopt in fully 
GR simulations of binary neutron star mergers. In §3, we summarize the methods 
used for analysis of gravitational waves. In §4, numerical results arc presented, 
paying particular attention to gravitational waveforms. Section 5 is devoted to a 
summary. Throughout this paper, we adopt geometrical units in which G = c = 1, 
where G and c are the gravitational constant and the speed of light. Latin and 
Greek indices denote spatial components (1-3) and space-time components (0-3), 
respectively. Sij{= 5^^) denotes the Kronecker delta. We use Cartesian coordinates 
x'' = {x, y, z) as the spatial coordinates and t = x^ denotes the time coordinate. 

§2. Formulation 

2.1. Basic equations 

We solve the Einstein and GR hydrodynamic equations without any approxi- 
mation by numerical simulation. Our formulation for a numerical solution of these 
coupled equations is described in detail in Refs. 18), 19) and 16), and therefore, we 
here only briefly review the basic equations. 

We write the line element in the form 

ds^ = g^ydx^dx" = {-a^ + PkP'')dt'^ + 2Pidx'dt + -/ijdx'dx^, (2-1) 

where gf^^y, a, f3^ (/3j = Jij^-'), and jij are the 4D metric, lapse function, shift vector, 
and 3D spatial metric, respectively. Following Refs. 18), 19) and 16), we define the 
quantities as 

7 = det(7i,) = e^^t (2-2) 
jij = e-^%j, (2-3) 

Aij = e-^t'(^Kij-^jijKy (2-4) 

where Kij is the extrinsic curvature, and K its trace. With this definition, det(7ij) is 
required to be unity. In the numerical computations, we evolve 0, ^ij, K, and Aij in 
time, instead of jij and K^j. We note that the indices of A^j (A^^) are raised (lowered) 
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in terms of 7*-^ (lij)- Hereafter, we use Di and Di as the covariant derivatives with 
respect to jij and jij, respectively: We define A = D^Di and A = D^Di. 

As the matter source of the Einstein equation, we adopt a perfect fluid, for which 
the energy-momentum tensor is written 



Tun 



(2-5) 



where p, s, P, and are the baryon rest-mass density, specific internal energy 
density, pressure, and four-velocity, respectively. We give initial conditions using 
polytropic equations of state as 



Kp 



r = i + 



1 



n 



(2-6) 



where k and n are a polytropic constant and a polytropic index. In the numerical 
time evolution, we use a T-law equation of state of the form 



P={r- l)pe. 



(2-7) 



In the absence of shocks, the polytropic form of the equations of state is preserved. 
We set n = 0.8 {F = 2.25) and 1 (i"' = 2) as a reasonable qualitative approximation 
to moderately stiff equations of state for neutron stars. 

The hydrodynamic equations [continuity, Euler and energy (or entropy) equa- 
tions] are written in the forms 



dtp* +di{p*v') = 0, 
dt{p*Uk) + diip^Ukv') -- 

dte* + di{e*v^) = Cart, 



-ae 



■.{P + Pe^)-P* 



+- 



ae 



-^^'^ 



-dk. 



w 



where = d/dx^^ p* = pwe^'^, h = 1 + e + P/p, w = au* , 
{pef/^we^^, and 



(2-8) 

, (2-9) 
(2-10) 

(2-11) 



Here, Part and ^a.rt are artificial viscous terms. In numerical simulation, we solve 
Eqs. (2-8)-(2-10) to carry out the evolution of p*, Uk and e*. 

Once Ui is obtained, w is determined from the normalization relation of the 
four- velocity, u'^u^ = — 1, which can be written as 



w 



1 + e-^'^f^itiUj 



1 + 



Pel 



1 -2 



p*{we^'t>) 



r-i 



(2-12) 
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The Einstein equation is split into constraint and evolution equations. The 
Hamiltonian and momentum constraint equations are written as 



AijA'^ + -K^ = 167rpH, 



or, equivalently 



-DjK = SnJj, 



(2-13) 
(2-14) 

(2-15) 
(2-16) 

where = e'^, pa = T^^'^n^riy and Jj = —T^^n^'yyi with = (— a,0). Here, 
Rij denotes the Ricci tensor with respect to 7jj, and Rj^ the Ricci scalar. These 
constraint equations are solved only at t = to set initial conditions (sec §2.3). 

Following Refs. 18), 19) and 16), we write the evolution equations for the geo- 
metric variables as 



AiIj = ^Rf!' - 27rpH^^ - ^ AijA'^ - -K 



{dt - P^di)jij = -2aAij + ^ikp'^j + ^jk/S^^i - -%(5^^k, 



(2-17) 



{dt-f3'di)A 



V 



-40 



ai - le^^iiRj" 



HJ 3 



k\ , ok 



1 



2 



+aiKAij - 2AikA/) + P^^^Akj + P^jAki - ^P^^kAj 



idt-l3'di)K = a 



- Aa + 47ra(pH + 5"/), 



(2-18) 
(2-19) 
(2-20) 



denotes the partial derivative and Sij = T^^^^^i'^yj. We solve these 



where " , 

evolution equations to carry out the evolution of 7jj, A^j, cp and K. 
For the computation of Rij, we decompose 

Rij = Rij + Rfj ) 

where Rij is the Ricci tensor with respect to jij, and 

Rfj = -2DiDjcj) - 2jijA(j) + 4Di<pDj<j) - 4jijDk<PD''q 
Rij is written as 



R 



5^''{—hij^kl + hik,ij + hjk,ii) + 2dk{f'"'ri^ij) - 2rLrj^ 



kl ; 



(2-21) 



(2-22) 



(2-23) 
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where we split "fij and 7*-' as 6ij + hij and 6^^ + /'•' , respectively. F^j is the Christoffel 
symbol with respect to ^ij, and r^^ij = %irlj. Because of the definition det(7ij) = 1, 
we have = 0. 

In addition to the Laplacian of hi j, R-ij involves terms linear in hij as S'^^'hik^ij 
and S^^'hjk^ii. To perform numerical simulation stably, we introduce an auxiliary 
variable Fi = S^^dijij,^^^ which evolves according to the evolution equation 

{dt - (3'di)Fi = -167raJi + 2al^f''^Ak,j + f'jAik - \A^%i,i + Q<t>,kA\ - \k,^ 

+ 5^^i^-2akA,^+l5\kh,,,i + {iu^\,+i,il3\, - J-(2.24) 

In the numerical simulations, we evaluate d^^'^ik^ij in Eq. (2-23) as Fi^. 

In the evolution equation for Aij, we compute Rj^ directly as RijY^ = e^'^'^ [Rij -j- 
Rfj)j''^ to suppress violation of the condition Jl/ = 0. We do not substitute the 
Hamiltonian constraint equation in this case because the Hamiltonian constraint is 
not exactly satisfied in numerical simulations, and its violation could systematically 
accumulate to eventually violate the condition A/ = significantly. 

In addition to the above formulation for the basic equations, we adopt a trans- 
formation for computation of Rj^^. In computing Rj^'' = 'y^^Rij, there exist two 
linear terms in hij, S'^^'hij^kl and j^^Fij. Although they are linear in hij, these 
terms should be small outside the strong field zone, because they are composed of 
the trace part of hij (S^^hij) and a vector component of hij (Fi), respectively. The 
magnitude of the term associated with Fi is controlled by the choice of the spatial 
gauge condition, and in our choice of the gauge (see §3), it is maintained to be small 
outside the strong field zone. As shown in (2-25), the other term is essentially the 
nonlinear term and should be small. However, in numerical calculations, the non- 
linearity is not guaranteed in the form 'y^^S'^^'hij^^i. Thus, to guarantee the nonlinear 
nature of this term, we rewrite it as 

rS^'hij^ki = -S'''hij,kf\i. (2-25) 

As a result of these procedures, Rf!' is maintained to be small outside the strong 
field zone. 

2.2. Gauge conditions 

We adopt the approximate maximal slice (AMS) condition if pa as the time 
slicing condition, and an approximate minimal distortion (AMD) gauge condition as 
the spatial gauge condition. 

The equation imposing maximal slicing, K = = dfK, is set up from the right- 
hand side of Eq. (2-20) as 

Aa = a (^AijA'^ + ^K^"^ + 47ra(pH + S^'') = 47rSa. (2-26) 

In Eq. (2-26), we leave K on the right-hand side, because it is not possible to 
guarantee it to be exactly zero in numerical simulations. 
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To impose the AMS condition, we solve the following parabolic-type equation 
for In a at each timestep, until an approximate convergence is achieved: 

dr \na=A\na+ (Djt In a) {D^ In a) - 47r(pH + S^^) - AijA'^ - + faKpl^^ . 

(2-27) 

Here r is a control parameter, and fa is a constant for which we take to be of 0(1). 
In solving Eq. (2-27) at each timestep, we substitute a trial function of a, which 
is extrapolated as Q!(x*) = 2a_i(x*) — a_2(a;*)) where a_i and a_2 sue the lapse 
functions at previous timesteps, and then carry out iteration to achieve convergence. 
We typically carry out this iteration for about 30 steps. 

The last term of Eq. (2-27) is introduced as a driver to achieve K = 0: Assuming 
that convergence is achieved and that the right-hand side of Eq. (2-27) becomes zero, 
the evolution equation for K can be written 

{dt - P^di)K = -j^aKpT. (2-28) 

Thus, if K is zero initially and convergence is completely achieved, the maximal 

slicing condition = is preserved. Even when the convergence is incomplete and 
K deviates from zero, the right-hand side of Eq. (2-28) causes \K\ to approach to 

— 1/2 ~i 

zero on the local dynamical timescale ~ .*-' Hence, the condition X = is 
expected to be satisfied approximately. 

To impose the AMD gauge condition, we solve the simple elliptic-type equations 

^flat^'i = Si, Z\flatr? = -Six\ (2-29) 

where /^flat is the Laplacian in the flat 3D space, and 

Si = 167ra Ji + lAij {&a - 6a&(f)) + ^aDiK. (2-30) 

The equations for Pi and rj are solved under the outer boundary conditions as 

where Cxx,Cxy,Cyx,Cyy,Czz, and are constants computed from the volume in- 
tegrations of SiX^ and S'jX*. Here, to derive these boundary conditions, we assume 
that the system is composed of two identical neutron stars. 
Prom Pi and rj, we determine /?' as 



lPi-liri,i + Pk,ix'') 



(2-31) 



*^ In the last term of Eq. (2-27), another function in place of p*, e.g., So,, might be a better 
choice as a driver to achieve K = 0. Other methods may be tried in the future. 
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Hence, /3* satisfies an elliptic-type equation of the form 

1 
3 



SijAR,tf3' + ^f3\j = Sj. (2-32) 



As described in Ref. 19), if the action 

1 = 1 d'x{daij){dt%i)f'j^' (2-33) 

is minimized with respect to we obtain the equation of a minimal distortion gauge 
condition ^^)*) for as 

7,fe4/3^ + ^DjOiP' + RjkP' = Sj. (2-34) 

Thus, the equation for in the AMD gauge condition is obtained by ignoring the 
coupling terms between and hij in Eq. (2-34). Since the ignored terms arc 
expected to be small, we can expect that / is approximately minimized in the 
AMD gauge condition. By ignoring these coupling terms, the equation for /3* is 
significantly simplified, enabling us to impose the spatial gauge condition with a 
relatively cheap computational cost. 

The other benefit of the AMD gauge condition is that Fi is guaranteed to be 
small everywhere except in the strong field region just around a highly relativistic 
object. This implies that the transverse condition S'^Wijjk = approximately 
holds for in the wave zone, helping extraction of gravitational waves near the 
outer boundaries of the computational domain. 

One drawback of the AMD gauge condition together with the maximum slicing 
is that the resolution quickly becomes bad in high density regions whenever the 
matter source collapses to a black hole. To improve the resolution, we modify the 
AMD gauge condition around a black hole formation region, subtracting the radial 
part of 13'' with a method described in Refs. 19), 13) and 14). 

2.3. Solution for initial conditions 

Even just before the merger, binary neutron stars that are not extremely compact 
are in a quasiequilibrium state because the timescale of the gravitational radiation 
reaction at Newtonian order ~ 5/{64J7(MNi?)^/^} (where M-^ and f2 denote the 
Newtonian total mass of system and the orbital angular velocity of binary neutron 
stars) is several times longer than the orbital period. Thus, for a realistic simulation 
of the merger, we should prepare a quasiequilibrium state as the initial condition. 
In simulations carried out to this time, we have constructed such initial data sets in 
the following manner. 

First, we assume the existence of a helical Killing vector, 



IT.] +^ hrJ ■ (2-35) 




*^ Definition of the minimal distortion gauge in this paper is shghtly different from the original 
version in Ref. 21). 
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Since the emission of gravitational waves violates the helical symmetry, this assump- 
tion does not strictly hold in reality. However, as mentioned above, the emission 
timcscalc of gravitational waves is several times longer than the orbital period, even 
just before the merger (c/. Table I), so that this assumption is acceptable for com- 
puting an approximate quasiequilibrium state. 

In this paper, we consider irrotational binary neutron stars with polytropic equa- 
tions of state. The assumption of an irrotational velocity field with the existence of 
the helical Killing vector yields the first integral of the Euler equation, ^"^^ 

h , ^^i. 
—7 + hukV = const, 
It* 

where = v'' — i'' . Because of the irrotational velocity field, Ui can be written as 

Ui = h~^Di^, (2-37) 
where # denotes the scalar velocity potential that satisfies the elliptic-type PDE 

A ipah-^D'<P) - Di [pah-\e + = 0. (2-38) 

The details of numerical methods for obtaining <P are given in Ref. 25). 

Initial conditions for geometric variables are obtained by solving the constraint 
equations (2-13) and (2-14) and equations for gauge conditions. Currently, we restrict 
our attention only to initial conditions in which hij = = dthij and K = 0. 

Using the conformal factor ^ = e"^, Aij = ip^Aij and A^^ = ip^A^^ , the Hamilto- 
nian and momentum constraint equations can be rewritten in the forms 

^flatV' = -SttphV-' - IAjA'^^-^ = S^, (2-39) 

A.\- = SnJi'iP'^. (2-40) 

The momentum constraint equation is further rewritten as elliptic- type equations 
using two methods. In one method, we use the York decomposition^^) as 

Aij = Wij + Wj,i - ^SijS^'Wk,i. (2-41) 

Then, denoting Wj as 

W^ = lBi-lix,i + Bk,^x''), (2-42) 
o o 

where x ^re auxiliary functions, Eq. (2-40) can be decomposed into two 

simple elliptic-type equations, 

Afi^tBi = SirJii^^ zdflatx = -SnJix'^^. (2-43) 

Since Jitp^ (= p>^Ui) is nonzero only in the strong field region, the solution of the 
momentum constraint equation can be accurately obtained even when outer bound- 
ary conditions are imposed near the strong field zone. This method is in particular 
useful in the case that Jitp^ is a priori given. 



(2-36) 
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In the second method, we write Aij in terms of /3* as 

Aj = ^ (sikP", + - \S^13\^ , (2-44) 

where we use hij = = dfhij aX t = 0. Then, the equation for is 

6ijA{i^tP' + ldidjP^+djln(^^^ (^diP^+6ikd^'dip''-'^6{dkP'^ = IGiraJi. (2-45) 

Thus, Aij can be obtained from either Eq. (2-41) or Eq. (2-44). Quasiequihbrium 
states are obtained by solving equations for instead of those for Wi, because 

we do not know Jitp^ a priori in this case, and also because we need to obtain /?* to 
solve the hydrostatic equations. Thus, to obtain a quasiequihbrium state, we solve 
Eqs. (2-36), (2-38), (2-39), (2-45), and (2-46) (see below) iteratively until sufficient 
convergence is achieved. 

In the following, modifying Uj for a fixed density configuration, we slightly re- 
duce the angular momentum and/or add an approaching velocity to quasiequihbrium 
states at t = to slightly accelerate the merger. Whenever we add such perturba- 
tions, we recompute the equations for Wi and ijj for the perturbed values of Ui to 
guarantee that the constraint equations are satisfied at t = 0. Then, to adjust 
the gauge conditions, we also recompute Eq. (2-45) and then Eq. (2-46) for given 
p*,Ui, Kij and ■0. 

In addition to the constraint equations, we solve an elliptic-type equation for a 
to impose K = = dfK initially. In the conformally flat 3D space, this equation is 
written 

^flat(aV') = 27raV''(pH + 23^'') + ^a^p-^AijA'^ = S^^. (2-46) 
Equations for ip, aip, and x ^I'S solved under outer boundary conditions as 

V' = l + ^ + 0(r-3), aV' = l + ^ + 0(r-3), = ^^^^fi^ + 0{r-% 

where C^, C^^, C^^,, C'^y, Cy^, C'yy, C'^^, and C-^ are constants that can be computed 
from the volume integrations of S^, Satp, JiX^ and JiX^. Here, to derive these bound- 
ary conditions, we assume that the system is composed of two identical neutron stars. 

2.4. Definitions of quantities 

In numerical simulations, we often refer to the total baryon rest-mass, ADM 
mass and angular momentum of the system, which are given by 

M^= J d^xp^, (2-47) 
Madm= f D'l/jdSi 
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/ 



J= 



1 

Stt 



'Aje^'^'dS, 



I 



r— >oo 

d^xe^'^ 



1 



(2-48) 



(2-49) 



where dSj = r'^Djrd{cos6)dip and ip^ = —y{dxy + x{dyy . Here, is a conserved 
quantity, and it uniquely specifies a model of stable binary neutron stars for a given 
r. Madm and J decrease as a result of the radiation reaction of gravitational waves. 

At t = 0, the ADM mass and angular momentum are given in the confer mal fiat 
three space with K = as 



-^ADMO 



d^x pnip^ + 



167r^ 



Jo = I d^xJi^'i}^. 



(2-50) 
(2-51) 



From these quantities, we define a nondimensional angular momentum parameter as 
q = Jo/-/Wadmo- 

Instead of , we specify a model of binary neutron stars using the compactness 
(M/i?) oo, which is defined as the ratio of the ADM mass to the circumferential radius 
of a spherical neutron star in isolation (see Tables I and II). To specify a model, we 
also use the ratio of the total baryon rest-mass of the system to the maximum allowed 
mass of spherical neutron stars for a given equation of state M*max, 



ih ■ 

*max 



(2-52) 



In addition to the above quantities, we use the relation between the baryon 
rest-mass and specific angular momentum that we define as 



M,(j) = / d'^x'p.ix'), 
Jj'>j 



(2-53) 



where J denotes the specific angular momentum j = hui(p^. Here, M^,{j)/M^ denotes 
a baryon rest-mass fraction whose specific angular momentum is larger than a value 
j (i.e., M,(j = 0)/M, = 1). 

Physical units enter the problem only through the polytropic constant k, which 
can be chosen arbitrarily or else completely scaled out of the problem. Since k"/^ (in 
geometrical units c = G = 1) has the dimension of length, dimensionless variables 
can be constructed as 



M,K-"/2, Madm = Madmk""/', R = Rk-''/'', 



J = Jk 



and Q = Qk^I"^ . 



(2-54) 



In the following, we present only these dimensionless quantities. 
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2.5. Code tests and accuracy estimation of simulations 

Several test simulations, including the spherical collapse of dust, stability of 
spherical stars, mode analysis of spherical stars, and evolution of rotating stars, 
have been carried out to check the accuracy of the numerical code. A list of these 
test simulations and their results are given in Ref. 16). (See also Ref. 22).) 

For the simulations presented in this paper, we monitored the violation of the 
Hamiltonian constraint, conservation of the baryon rest-mass, conservation of the 
ADM mass [Madm+ (energy loss through gravitational radiation) should be con- 
served] and conservation of angular momentum [J+ (angular momentum loss through 
gravitational radiation]. Simulations were stopped when the conservation of the 
ADM mass was significantly violated. 

As described in Ref. 13), the violation of the Hamiltonian constraint is locally 
measured by the equation as 



u = 



(2-55) 



We briefly discuss the magnitude of and conservation of the mass and angular 
momentum in the Appendix. 

§3. Analysis of gravitational waves 

As in Ref. 19), we analyze gravitational waveforms using two methods. In both 
cases, we extract gravitational waves near the outer boundaries of the computational 
domain. 

In one method, we compute the following quantities along the z-axis: 

^+ ^ 2Mad Jo7m/^)oo " ^ MADMo(M/i?)oo'^^^- ^^'^^ 

Here ^obs denotes the location at which we extract gravitational waves. Since the 
gauge condition we adopt is approximately transverse and traceless in the wave zone, 

and hx arc expected to be appropriate measures of gravitational waves emitted. 
The amplitude of gravitational waves will be largest along the z-axis as a result of 
our choice of the orbital plane, so that h+^x can be used to measure the maximum 
amplitude of gravitational waves. We also note that and are composed only 
of jmj = 2 modes, because other modes vanish along the z-axis. 

The amplitude of gravitational waves for equal-mass binaries along the z-axis in 
the quadrupole formula is 

^GW = — (3-2) 

where Mn is a Newtonian mass and a an orbital separation. At the point of contact 
of two spherical stars in a circular orbit, we have 

ftow = (3.3) 

■^-obs 
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Thus, by normalizing the amphtude by Madmo and (M/R)^, as in Eq. (3-1), the 
maximum values of h+ and are scaled to become approximately unity. 

To search for the dominant frequency of gravitational waves, we compute the 
Fourier spectra defined as 



h+,xif) 



(3-4) 



In the analysis, we choose to be the time at which wc stop simulation, for simplicity. 
Before t < Zohs, no waves propagate to Zobs, so that we choose t'l ~ ^ohs • 

Gravitational waves are also measured in terms of the gauge invariant Moncrief 
variables in a flat spacetime in the following manner. ^^)*) First, we carry out a co- 
ordinate transformation for the three-metric from Cartesian coordinates to spherical 
polar coordinates, and then split as r]ij + X^im Cij") where rjij is the flat metric in 



spherical polar coordinates and C,\^ is given by 



Aim 

5jj 



/ H2l'mYlm hiimYi 

rriM 



hllmYlm,<p 

r^sm^e{KimYi^-GimWim) J 



+ 



/ 



V * 



-Cimd^Yim/ sin 9 
r^Di^Xi^/ sine 
* 



CimdeYim sin 9 
-r'^Di^Wim sin 9 
-r^DimXim sin 9 



(3-5) 



Here, * denotes the relation of symmetry. The quantities i?2im) hum, Kij^, Gim, 
Cijn and Di^ are functions of r and t, and are calculated by performing integrations 
over a two-sphere of given radius [see Rcf. 18) for details]. Also, Yj^ is the spherical 
harmonic function, and Wi^ and Xi^^ are defined as 



{def - cot 9d0 - 



1 



sm 



Y. 



Imi 



dg — cot 9 



Yim- (3-6) 



The gauge invariant variables of even and odd parities are then defined as 



where 



kllm= Kim + 1(1 + l)Glm + '^rdrGim - 2 



llm 



k2lm= 



H2lm _ 1 d 

2 2dr 



r{Kim + + l)Gim} 



(3-7) 
(3-8) 

(3-9) 
(3-10) 



*^ The Moncrief formalism was originally derived for the Schwarzschild spacetime. We here 
apply his formalism in a flat spacetime. 
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Using the above variables, the energy luminosity and angular momentum flux 
of gravitational waves can be calculated as 



dE _ ^ 

l,m 

dJ _ 

l.m 



\midtRfjRfj + \m{dtR?jR?J 



(3-11) 
(3-12) 



The total radiated energy and angular momentum are calculated from 

dE _ r dJ 



AE 



Wc have computed modes with / = 2, 3 and 4, and found that even modes with 
I = \m\ = 2 are dominant. For this reason, in the following, we focus only on these 
modes. 

§4. Numerical results 

4.1. Products after merger 

In Table I, we list several quantities that characterize qTiasicquilibrium states of 
irrotational binary neutron stars used as initial conditions in the present simulations. 
All quantities are dimensionless and appropriately scaled with respect to k. We 
choose binaries at innermost orbits for which the Lagrange points appear at the inner 
edge of neutron stars. Since the orbits of such binaries are stable, we reduce the 
angular momentum slightly to accelerate the merger. Also, we add an approaching 
velocity at t = for some cases (see discussion below). The values of Madm and 
q listed in Table I are those of quasiequilibrium states before reducing the angular 
momentum and before adding the approaching velocity. Note that the frequency of 
gravitational waves for these quasiequilibria are given by 

f ^0 ^.n.nx. f 2-8Me V Co 

/QE = -- 960 Hz(^^^j(^^j , (4-1) 

where f2o denotes the angular velocity of quasiequilibrium states. Thus, the orbital 
period of the quasiequilibria, Pt=o, is 



/ 2.8M0 \ V Co ^ 



(4-2) 



Here, Cq is a compactness parameter of orbits defined as 

Co ^ (Madmo^^o)^/^ - (4-3) 



where a denotes the initial orbital separation. Products of the merger which we 
found when we stopped simulations are also described in Table I, in the last column. 
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Table I. A list of several quantities for quasiequilibria of irrotational binary neutron stars with 
r = 2.25 and 2. The compactness of each star in isolation, {M/R)oo, the maximum den- 
sity, Pmax = ^"pmax, the total baryou rest-mass, M* = k~"^^M,, the ADM mass at i = 0, 
Madmo = k~"''^Madmo, q = Jo/M^omo, Pt=o = Pt=o/MAUMo, the orbital compactness, 
Co = (Madmo^^o)^''^, the ratio of the emission timcscale of gravitational waves to the orbital 
period at Newtonian order, Rj^ = 5(MADM0^2o)~°/^/1287r, the ratio of the baryon rest-mass of 
each star to the maximum allowed mass for a spherical star, Rmasa = Aft/M^^'max, and products 
we found when we stopped simulations are listed. All quantities are normalized by k appropri- 
ately to be dimensionless: We can rescale the mass to a desirable value by appropriately choosing 
K. Here, Af^^'max denotes the maximum allowed mass of a spherical star (M^^'^ax ~ 0.162 at 
Pmax « 0.524 for r = 2.25 and M^max « 0.180 at p^ax « 0.32 for T = 2). BH and NS mean 
"black hole" and "neutron star". MA (marginal) implies that we were not able to determine 
the product. 



Model 


r 


(M/R)^ 


Pmax 




Madm 


q 


Pt=o 


Co 


Rr 


-^mass 


Product 


(A) 


2.25 


0.12 


0.139 


0.186 


0.173 


1.02 


233 


0.0897 


5.1 


1.15 


NS 


(B) 


2.25 


0.14 


0.169 


0.216 


0.198 


0.97 


182 


0.106 


3.4 


1.33 


BH 


(C) 


2.25 


0.16 


0.202 


0.244 


0.220 


0.93 


145 


0.124 


2.3 


1.51 


BH 


(D) 


2.25 


0.17 


0.220 


0.257 


0.231 


0.92 


130 


0.132 


2.0 


1.59 


BH 


(E) 


2 


0.10 


0.0695 


0.224 


0.211 


1.07 


315 


0.0735 


8.5 


1.24 


NS 


(F) 


2 


0.11 


0.0798 


0.242 


0.227 


1.03 


271 


0.0813 


6.6 


1.35 


MA 


(G) 


2 


0.12 


0.0910 


0.261 


0.243 


1.00 


236 


0.0892 


5.2 


1.44 


BH 


(H) 


2 


0.14 


0.117 


0.292 


0.270 


0.94 


184 


0.105 


3.4 


1.62 


BH 


(I) 


2 


0.16 


0.149 


0.320 


0.292 


0.90 


146 


0.123 


2.3 


1.78 


BH 



4.1.1. Set-up for computational grids 

We have performed the simulations using a fixed uniform grid and assuming 
reflection symmetry with respect to the z = plane, (i.e., the equatorial plane is 
chosen as the orbital plane.) The typical grid size is (505, 505, 253) for {x,y,z). 
The largest grid size is (633, 633, 317). To investigate numerical efi'ccts associated 
with the location of outer boundaries, we performed simulations choosing smaller 
grid sizes. The grid covers the region —L<x,y<L and < z < L, where L 
is a constant. The grid spacing is determined from the condition that the major 
diameter of each star is covered with 33 grid points initially. To investigate effects 
of numerical dissipation and diffusion, we also choose larger grid spacing in the 
calibration. With a (505, 505, 253) grid size, the computational memory required is 
about 120 GBytes, and the computational time for one model is typically about 100 
CPU hours for about 10000 timesteps using 32 processors on FACOM VPP5000 in 
the data processing center of NAOJ. 

To perform an accurate simulation excluding spurious effects associated with 
imposing outer boundary conditions in a local zone, L should be larger than the 
wavelength of gravitational waves, A. Using almost the entire capacity of the super- 
computer VPP-5000 at NAOJ, it is possible to extend L as large as ^ A. However, 
such a large simulation is limited by the computational time assigned to us. Since the 
grid spacing should be smaller than 1 /30 of the diameter of each star to avoid large 
numerical dissipation and diffusion (see below) , the size of the computational domain 
along each axis, L, becomes about 0.5-0.6 Aq in this work. (Here, Aq = vr/l2o-) This 
implies that gravitational waves in the early stage are not very accurately computed. 
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Table II. List of the various parameter values in the simulation for models (A), (B), (C), (D), (E), 
(F), (G), (H) and (I): The model, the compactness, the q parameter at t = 0, the ratio of the 
approaching velocity to orbital velocity, the grid number, L in units of Aq, 1/ in units of Madmo, 
and the grid number inside the initial gravitational mass length Madmo/^s; are listed. Here, 

v„b = (MadmoJ^o)'/^ = C(y^ 



A/Tnrlpl 


(M IU\ 


Q 


^app / ^orb 


1 -.VI (H n n TY^ r">0"r 




L/ A/adm 




/ A 

(A-UJ 


n 1 o 
U. IZ 


i.UU 


A 
U 




n Qoo 
U.ozz 


0.7 ^ 


Q on 
o.yu 


( ^ ^\ 


n 1 o 
U. iz 


1 nn 
I.UU 


n 1 
U.i 


(Oio, OiO, ID i ) 


n Q/1 /I 




o.yu 


[A-Z) 


n 1 o 
U. iz 


i.UU 


n 
U 


l^oUo, oUo, Zoo) 


n KKR 

U.ooO 


dA 7 


Q on 
o.yu 


(A-.j) 




1 nn 

i.UU 


u 




U.OOO 


f\A 7 


9 01 
z.y 1 


(B-0) 


0.14 


0.95 





(293, 293, 147) 


0.343 


31.1 


4.69 


(B-1) 


0.14 


0.96 


0.1 


(313, 313, 157) 


0.366 


31.1 


4.69 


(B-2) 


0.14 


0.95 





(505, 505, 253) 


0.592 


53.7 


4.69 


(B-3) 


0.14 


0.95 





(377, 377, 189) 


0.592 


53.7 


3.50 


(C-0) 


0.16 


0.91 





(293, 293, 147) 


0.364 


26.3 


5.55 


(C-2) 


0.16 


0.91 





(505, 505, 253) 


0.628 


45.4 


5.55 


(C-3) 


0.16 


0.91 





(633, 633, 317) 


0.788 


56.9 


5.55 


(D-0) 


0.17 


0.89 





(313, 313, 157) 


0.398 


26.0 


6.00 


(E-1) 


0.10 


1.06 


0.01 


(505, 505, 253) 


0.492 


77.6 


3.25 


(F-1) 


0.11 


1.02 


0.01 


(505, 505, 253) 


0.514 


69.5 


3.62 


(G-0) 


0.12 


0.99 


0.01 


(249, 249, 125) 


0.263 


30.9 


4.01 


(H-1) 


0.14 


0.93 


0.02 


(505, 505, 253) 


0.569 


52.3 


4.82 


(I-l) 


0.16 


0.89 


0.02 


(505, 505, 253) 


0.606 


44.1 


5.71 



On the other hand, the wavelength of quasi-periodic waves of the merged object ex- 
cited during merger is much smaller than Aq and L, so that the waveforms in the 
merger stage can be computed accurately. This point is discussed in detail in §4.2. 

4.1.2. Set-up for initial data 

As found in Ref. 5) , orbits for all irrotational binaries of equal mass with F < 
2.5 are dynaniically stable from infinite separation to the innermost orbit at which 
Lagrange points appear at the inner edges of neutron stars. Thiis, the merger in 
reality should be triggered by the radiation reaction of gravitational waves. As 
argued above, however, the radiation reaction in early stages during which A > L 
is not accurately computed with our current computational resources. To induce 
prompt merger by destabilizing the orbital motion, we slightly decrease the angular 
momentum at i = from the equilibrium value (by less than 1%).*^ 

*^ As initial data sets for the simulation with F = 2.25 and for the simulations in Ref. 13), we 
adopted solutions of binary neutron stars in quasiequilibria which arc computed using the numerical 
code used in Ref. 5). Our latest computation for the quasiequilibrium computed with a revised 
code indicates 2-3 % smaller angular momentum than in the previous computation, while other 
quantities, such as the gravitational mass and angular velocity, do not vary much. In our previous 
code, the magnitude of the shift vector is slightly overestimated around the strong field region. As a 
result, the angular momentum is systematically overestimated. We became aware of this fact after 
we finished most of the simulations for F = 2.25. Thus, all the results presented in this paper were 
computed using the old initial data. We have adopted the new data sets only for F — 2 cases. The 
relevant quantities for F = 2.25 computed using the new code are as follows (compare with Table 
I)- 
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Using the quadrupole formula, the angular momentum loss in one orbital period 
z^Jone is estimated as 47rMADMoC'o/^' ^"^^ hence the ratio of AJone to the total 
angular momentum J is estimated as 



J 




(4-4) 



Thus, the artificial decrease of the orbital angular momentum by 1% is much smaller 
than the loss through gravitational radiation in one orbital period. 

As wc mentioned in Rcf. 13), we also performed several test simulations with 
small grid number varying the depletion factor in the range 2%, but we found that 
the qualitative results for the merger stage are not modified at all. The depletion only 
affects the evolution of the late inspiraling stage before the hydrodynamic interaction 
of two neutron stars sets in. 

To investigate the effect of the approaching velocity on the product and gravi- 
tational waveforms, as well as to accelerate the merger in some cases, we also add 
an approaching velocity v^pp for some models by changing the component of the 
four-velocity as 



X (^a;) original 



^app 



(4-5) 



where we assume that the center of mass of two stars is located along the x-axis 
at t = 0. For F = 2, we add Va,pp as in Ref. 13). Explicitly, we set v^pp = 



1/2 

• 



(0.01 — 0.02)7;orb) where the orbital velocity forb is defined as (Madmo^o)^^'^ = C^ 
The order of magnitude of this approaching velocity is in approximate agreement 
(within a factor of 2-3) with that induced by the radiation reaction of gravitational 
waves. For F = 2.25, v^pp was basically chosen to be zero. To investigate its effect 
for comparison, we choose a large nonzero value for some models as Wapp = 0.1i;orb 
(see §4.2.3). 

These treatments of the initial conditions could cause a certain systematic de- 
viation from a realistic merger in the early stages of the simulation, i.e., before the 



Model 


{M/R)oo 


Pmax 


M, 


Madmo 


<? 


-Pt=()/AfADMO 


Co 


(A) 


0.12 


0.139 


0.186 


0.173 


1.00 


232 


0.0901 


(B) 


0.14 


0.169 


0.216 


0.198 


0.95 


181 


0.106 


(C) 


0.16 


0.202 


0.244 


0.220 


0.91 


145 


0.123 


(D) 


0.17 


0.221 


0.257 


0.231 


0.89 


131 


0.132 



The latest computation is expected to provide more accurate initial data sets. Hence, reader 

might think it desirable to re-perform all simulations with the new initial data for F = 2.25. However, 
we think it not necessary to do this at present. Our present main purpose is to demonstrate that it 
is possible to extract gravitational waveforms within ~ 10% error as well as to determine the fate 
of mergers and the dependence of this fate on the compactness of neutron stars and equations of 
state. For these purposes, slight inaccuracy in the initial data sets is not important. By performing 
several test sinmlations, we have confirmed that qualitative properties are not modified using the 
new data sets for F = 2.25 with a small depletion factor (0-2%) and small grid numbers. Even 
the quantitative results agree well with those presented in this paper for a particular choice of the 
depletion factor in this range. Therefore we can conclude that the qualitative features in the merger 
stage {t > Pt=o) depend very weakly on the depletion factor, as long as the magnitude is of order 
0.01. 
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Fig. 1. Snapshots of the density contours for p in the equatorial plane for model (A-2). The 

contour curves are drawn for p/po = 10""'^-', where po is the maximum value of po at t = 
(here po = 0.139/v~"), for j = 0, 1, 2, • ■ ■ , 10. The maximum density at each time step is found 
in Fig. 4. Vectors indicate the local velocity field {v^,v'"), and the scale is shown in the upper 
right-hand corner. Pt=o denotes the orbital period of the quasiequilibrium configuration given 
at t = 0. The length scale is shown in units of GMadmo/c^, where Madmo is the gravitational 
mass at t — 0. Note that t = 3.18Pt=o « 740Madmo- Each panel does not necessarily represent 
contours of a constant time interval. 



hydrodynamic interaction of the two neutron stars sets in. Obviously, a more realis- 
tic simulation taking into account the radiation reaction in early stages with a large 
computational domain or with improved wave extraction techniques''^) is desirable. 
However, as we see below, qualitative and even quantitative results for the merger 
stage are not affected strongly by these initial conditions, nor by the few percent 
artificial decrease of the angular momentum. 
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t=0 OOP t=0 51Ft t=l 02Pt„ 




-10 -5 5 10-10 -5 5 10-10 -5 5 10 



X / M,„„ X / X / M,„„ 

Fig. 2. The same as Fig. 1, but for model (B-2). The contour curves are drawn for p/po = W~°'^^ 
for J = 0, 1, 2, • • • , 10. Here, po = 0.52k~", which is the critical density of stability for spherical 
stars against gravitational collapse. The thick solid circle in the final panel indicates the apparent 
horizon. 

4.1.3. Characteristics of merger for F = 2.25 

In Figs. 1-3, we display the density contour curves and velocity vectors for p and 
v'^ at selected timesteps for simulations of models (A-2), (B-2) and (C-S).*-' Time 
evolutions of the maximum density for p and a at r = for models (A-2), (B-2), 
(C-3) and (D-0) arc also shown in Fig. 4. The maximum density and central value 
of the lapse function at timesteps selected in Figs. 1-3 are plotted in Fig. 4. 

In the merger for model (A), the product we found when we stopped the simu- 

*^ Animations for these simulations can be seen at 
http://www.esa.c. u-tokyo.ac.jp/~shibata/anim.html. 
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lation at t ~ 3.3Pt=o, where Pt=o = 27r/I?o, is a massive neutron star. However, the 
merged object does not promptly settle down to an axisymmetric rotating neutron 
star. Soon after the onset of the merger, the merged object develops spiral arms 
(see fourth panel of Fig. 1), because it is rapidly rotating, and hence appears to be 
dynamically unstable with respect to formation of bar and spiral arms. The spi- 
ral arms do not spread widely outward, but wind about the central core (see fifth 
panel). After the winding, spiral arms develop again (see sixth panel), and subse- 
quently wind about the core after a short time (sec seventh panel). After repeating 
this process several times, the merged object gradually settles into an ellipsoidal 
rotating star (see eighth and ninth panels). During these oscillations, gravitational 
waves of fairly large amplitude are excited over a long timescale. This gravitational 
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0.8 
S 0-6 

A 0-4 
^ 0.2 


12 3 

t / Pt=o 

Fig. 4. Evolution of the maximum value of p and a at r = for model (A-2) (dashed curves), 
(B-2) (long-dashed curves), (C-3) (solid curves) and (D-0) (dotted-dashed curves). Solid circles 
denote the time slices which we choose to generate the contour plots displayed in Figs. 1-3. 
Note that the merger starts at f ~ Pt=o for each model and that apparent horizons are formed 
when a{r = 0) becomes ~ 0.04-0.05. 

radiation eventually triggers the collapse of this massive neutron star to a black hole 
(see discussion in §4.2). 

The total baryon rest-mass of the binary for model (A) is ~ 16% larger than the 
maximum allowed value of a spherical star in isolation. Even with such a large mass, 
the merged object does not collapse to a black hole within a couple of dynamical 
timescales. As indicated in Fig. 1, the merger proceeds very gradually, because the 
approaching velocity at the point of the contact of two stars is not very large. Con- 
sequently, the shock heating does not appear to play an important role in merging. 
Indeed, the function P/p^, which represents the change of the entropy, increases by 
at most 10% from the initial value (i.e., k) around the central region of the merged 
object. This implies that the rotational centrifugal force plays an important role for 
supporting the self-gravity of such a massive neutron star. To illustrate the impor- 
tance of the rotation in supporting the large mass, the angular velocity along x and 
y-axes, defined as l^^/xl and /y\, are displayed in Fig. 5. These plots show that 
the formed massive neutron star is differentially rotating, and the magnitude of the 
angular velocity is of order of the Kepler velocity, i.e., !7Madm(-^^/Madm)^^^ = 0(1), 
where R is the characteristic radius of the merged object (?s 5Madm)- 

For models (B), (C) and (D), in which {M/R)^ > 0.14, a black hole is formed in 
the present simulations. For all the cases, we were able to locate apparent horizons 
in the last stage of the simulations using our apparent horizon finder. However, 
the formation process is slightly different for each of these models (c/. Fig. 4). For 
model (D), a black hole is formed soon after the first contact of the two neutron 
stars at t Pt=o (see Fig. 4). For model (C), the formation timescale is slightly 
longer. In this case, the merged object first experiences a bounce (see fifth and sixth 
panels in Fig. 3, and Fig. 4), and then collapses into a black hole. For model (B), the 
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Fig. 5. The angular velocity in units of -M^q,^jq along x and y axes at t/Pt=o = 3.18 for model 
(A-2). 




1 2 




12 3 



t / Pt=o 

Fig. 6. Evolution of the maximum values of p and a at r = for models (A-2) (solid curves), (A-3) 
(dashed curves), (B-2) (soUd curves) and (B-3) (dashed curves). 

formation timescale to a black hole is even longer. In this case, the merged object 
quasi-radially oscillates for several times after the first contact. Indeed, the lapse 
function at r = does not monotonically approach to zero, as shown in Fig. 4. The 
collapse toward a black hole seems to set in after the angular momentum is dissipated 
through gravitational radiation. (Indeed, the angular momentum decreases by a large 
amount of 0(10%) from the late inspiral to the merger stages; see §4.2 and Table 
III.) As we show in §4.2, the difference in the formation process of black holes is 
significantly reflected by the waveform of gravitational waves (c/. Fig. 9) and the 
Fourier spectra (c/. Fig. 13). 

To investigate the effect of resolution on numerical results, we compare the time 
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evolution of the maximum density and the central lapse function for models (A-2) 
and (A-3) and for models (B-2) and (B-3) in Fig. 6. In these comparisons, we fix 
the location of the outer boundaries, while the grid spacing for (A-3) and (B-3) is 
1.34 times larger than that for (A-2) and (B-2), respectively (see Table II). Figure 6 
indicates that convergence is achieved fairly well. It also shows that with lower 
resolution, the numerical dissipation of the angular momentum is larger, so that 
the merger happens earlier. Also, high density peaks are captured less accurately, 
because of larger numerical diffusion. This makes the formation timescale of a black 
hole longer. It is interesting to note that, if we use much lower resolution, model (B) 
might not result in black hole formation in several rotational periods of the merger 
object, because of a large numerical diffusion. This result indicates that to determine 
the criterion of black hole formation and the formation timescale accurately, we need 
to perform better-resolved numerical simulations. 

A noteworthy result for models (B), (C) and (D) is that the baryon rest-mass 
fraction outside the apparent horizon at its formation is less than 1% of the total 
rest-mass. This implies that the fraction of the baryon rest-mass in the disk around 
the black hole is very small. We expect that the mass fraction of the disk is much less 
than 1% of the total baryon rest-mass in these cases. In the following, we describe 
the reason for these results. 

In Fig. 7(a), we plot M^{j)/M^ as a function of j'/Madmo for models (A), (B) 
and (C) at t = 0. It is found that there is no fluid element for which j/Madmo > 1-6 
for any model. (We also plot this curve for i"' = 2 in Fig. 7(b). This indicates that 
the spectrum of the specific angular momentum depends very weakly on F. Thus, 
the following conclusion holds irrespective of the value of F and the compactness 
of the neutron stars.) This small value of the specific angular momentum of all the 
fluid elements at i = plays an important role for determining the final outcome. 

As we found in the present simulations [as well as in Refs. 13) and 14)] a 
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Fig. 8. Evolution of the maximum values of p and q at r = for models (E-1) (solid curve), (F-1) 
(solid curve), (G-0) (dashed curve), (H-1) (solid curve) and (I-l) (solid curve). Note that the 
merger starts at < ^ Pt=o for each model and that apparent horizons are formed when a(r = 0) 
becomes ~ 0.04. 



quite large fraction of the fluid elements are eventually swallowed into the black hole 
whenever it is formed [i.e., for (B)-(D) as well as (G)-(I) (see Table I)]. Gravitational 
radiation carries energy from the system, but it is likely ^ 1% of Madmo (see §4.2 
and Table III). Thus, the ADM mass of the formed black holes is approximately 
equal to the initial value of the system. In contrast to the much smaller value for the 
energy, the angular momentum is dissipated by gravitational waves by ~ 10% of the 
initial value before the formation of black holes (see Table III). These facts imply that 
q = J/Mp^Y)^ of formed black holes is not equal to the initial value, ~ 0.9-0.95 but 
to ~ 0.8-0.85 for models (B)-(D). The specific angular momentum of a test particle 
in the innermost stable circular orbit around a Kerr black hole of mass Madm and 
q = 0.8 (0.9) is « 2.4Madm (2.1Madm)- Fluid elements with the specific angular 
momentum less than this value have to be swallowed into black holes. Therefore, 
Fig. 7 shows that no fluid element for irrotational binary neutron stars just before 
the merger has large specific angular momentum large enough to form a disk around 
the formed black hole. For a disk formation, a certain transport mechanism of the 
angular momentum, such as a hydrodynamic interaction, is necessary. Since the 
black holes are formed in the dynamical timescale of the system, this mechanism has 
to be very effective to increase the specific angular momentum of the fluid elements 
by more than 40 % in such short timescale. However, such a rapid process is unlikely 
to work, as indicated by the present simulations. 

4.1.4. Dependence of merger on F 

The results of the simulation for F = 2 are qualitatively the same as those 
for F = 2.25, shown in Figs. 1-6. As an example, we display the maximum value 
of p and a at r = as a function of time for models (E-1), (F-1), (G-0), (H-1) 
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and (I-l) in Fig. 8. Black holes are formed for models (G)-(I),*'* and a neutron 
star is formed for model (E). For model (F), a black hole is likely formed at some 
t ^ 3Pt=o ~ 800-Madmo- However, it is difficult to determine the time of formation 
accurately, because numerical error accumulates for t ^ 3Pt=o, resulting in inaccurate 
computation. We can be sure that at least the formation time would be ^ 3Pt=o, as 
in the case of model (A) . 

Comparison of Figs. 4 and 8 reveals that the time evolutions of the merger in 
models (E), (G), (H) and (I) for F = 2 are similar to those in models (A), (B), (C) 
and (D) for F = 2.25, respectively. A noteworthy point is that even if the value of 
illmass or {M/R) oo are identical, the evolution process of the merger depends on F; 
e.g., the compactness of the neutron stars in model (A) is identical to that in model 
(G). but the behavior of a at r = and Pmax differ for the two. The minimum value 
of R 

mass for prompt formation of a black hole (i.e., for formation of a black hole 
without any transient object) is ~ 1.5 for F = 2.25 and ~ 1.6 for F = 2. Thus, for 
softer equations of state (i.e., for smaller F), this value is larger. 

In terms of the compactness (M/R)^, the minimum value for prompt formation 
of a black hole is ^ 0.16 for F = 2.25 and ~ 0.14 for F = 2. Thus, for softer equations 
of state, this value is smaller. The compactness of a realistic neutron star of mass 
1.4M0 is likely in the range between 0.14 and 0.21, because the theory of neutron 
stars tells us that the radius is likely between 10 and 15km. "^^^ Thus, in a realistic 
merger, a black hole is likely formed soon after the onset of the merger for softer 
equations of state (with F ^2), while for stiffer equations of state (with F ^ 2.25), 
a transient object could be formed before forming a black hole in the merger of low 
mass binary neutron stars. 

4.2. Gravitational waveforms 

4.2.1. Character of waveforms and luminosity 

In Fig. 9, we plot and hx as functions of retarded time (t — 2obs)/-Pt=o for 
models (A-2), (B-2) and (C-2). (In the following, we always normalize the retarded 
time by Pt=o-) These results are obtained in a simulation with (505,505,253) grid 
number. The waveforms are quite similar in these three cases for times satisfying 
t — Zohs ^ Pt=o, but they vary as the merger proceeds toward the final outcomes. 

In the case of massive neutron star formation [model (A)], quasi-periodic gravi- 
tational waves of fairly large amplitude, which are excited due to nonaxisymmetric 
deformation and oscillation of the merged object, are emitted after the merger sets 
in. The wavelengths of these quasi-periodic waves are likely associated with several 
fundamental oscillation modes of the merged object, because the quasi-periodic os- 
cillation is not a simple sine-cosine curve, but has a modulation. Since the radiation 
reaction timescale is much longer than the dynamical (rotational) timescale of a 

In Ref. 13), wc concluded that the product after merger for model (G) is a massive neutron 
star. In the present work, we perform a simulation again for a longer timescale than the previous 
one, and find that a black hole is eventually formed as a result of collapse of a transient massive 
object. Even for the less massive models (E) and (F), a black hole will eventually be formed, due 
to dissipation of the angular momentum by gravitational radiation (see §4.2). Only the timescale 
for the formation of the black hole depends on the models used. 



27 




(t-z„J / Pt=o (t-z„J / Pt=o 



Fig. 9. h+ and hx as functions of the retarded time for models (A-2), (B-2) and (C-2). 



massive neutron star, the quasi-periodic waves will be emitted for many rotational 
cycles (sec below). 

Even for the case of black hole formation, quasi-periodic gravitational waves 
are excited due to the nonaxisymmetric oscillation of a transient merged object 
before collapsing to a black hole. Since the formation timescale of the black hole 
is different for models (B) and (C), depending on the initial compactness of the 
neutron stars, the duration of the emission of the quasi-periodic waves induced by 
the nonaxisymmetric oscillation of the merged objects is also different. Because the 
computation crashed soon after the formation of the apparent horizon, we cannot 
draw a definite conclusion from the present simulation with regard to gravitational 
waves in the last stage. However, we can expect that after the formation of black 
holes, QNMs of the black holes are excited and gravitational waves are eventually 
damped. To complete waveforms up to the ringing tail of QNMs, we need either 
to develop an implementation such as horizon excision or to use the perturbative 
extraction technique recently developed by Baker, Campanelli and Lousto. '^^^ 

The properties mentioned above are also reflected in the energy luminosity, as 
shown in Fig. 10. Here, we display the luminosity for I = \m\ = 2 modes as functions 
of the retarded time for (A-2), (B-2) and (C-2). We plot two luminosity curves that 
are extracted at robs ~ (solid curves) and 0.77L (dashed curves), respectively. It 
is found that during the early stage (i.e., for t— robs ^ Pt=o), the luminosity extracted 
at 0.77L is larger than that at L. This is because the wavelength of 

gravitational waves is longer than L during this stage, and hence gravitational waves 
are not accurately extracted from the metric. However, for t — robs ^ Pt=Oj the 
two curves agree fairly well, because the wavelength becomes shorter than robs > and 
consequently accurate extraction of gravitational waves becomes possible. We also 
note that it is difficult to obtain a smooth curve of dE/dt, because computing the 
energy luminosity requires taking time derivative of the gauge invariant quantity, 
which introduces a fairly large numerical noise. 
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Fig. 10. Luminosity of gravitational waves as functions of the retarded time for models (A-2), (B-2) 
and (C-2). Solid and dashed curves represent the luminosity computed at robs ~ L and 0.77 L, 
respectively, robs denotes the radius at which the gauge invariant variables are extracted. For 
i — fobs ^ Pt^o, the wavelength of gravitational waves is longer than L, but for t — robs ^ 1.2Pt=o, 
it becomes shorter than L (see Fig. 9). The unit of luminosity is c^/G ~ 3.63 x 10^^ erg/sec. 




We list the total radiated energy AE and angular momentum AJ in Table III. 
A typical magnitude of AJ/Jq is ~ 10%, which is consistent with the decrease of 
the angular momentum J computed by Eq. (2-49). (It was impossible to confirm 
the consistency of the change of Madm due to radiation reaction, because Madm 
varies by ^ 10% during simulations due to numerical error; see discussion in the 
Appendix.) It should be noted that AJ/Jq is much larger than AE/M^dmo, because 
the approximate relation A J AE / Q holds, and hence 



^lofiV ^ (4.6) 

{AE/Mabmo) \qj\MADMO^iJ' 



where J? here denotes a typical magnitude of the angular velocity. 

During the early stages (i ^ Pt=o), the luminosity of gravitational waves gradu- 
ally increases, reaching a peak. The peak luminosity is approximately proportional 
to (M/R)^, as expected from the quadrupole formula. However, after the peak is 
reached, the behavior of the curve depends strongly on the outcome of the merger. 



Table III. Total radiated energy and angular momentum in units of their initial values. We esti- 
mated the fluxes at robs ~ L. If we estimated them at robs ~ 0.771/, they are overestimated by 
~ 20%, because of inaccurate extraction of gravitational waves in the early stages during which 
the gravitational wavelength is longer than L. 



Model 


AE/Maumo 


AJ/Jo 


(A-2) 


0.5% 


8% 


(B-2) 


0.6% 


10% 


(C-2) 


0.9% 


12% 


(E-1) 


0.3% 


6% 


(F-1) 


0.3% 


7% 


(H-1) 


0.6% 


10% 
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For model (C-2), in which a black hole is formed in a short timescale at i ~ 2Pt=o 
(t — Tobs ~ 1.7Pt=o), a sharp peak is found before the black hole is formed. This peak 
is likely associated with an oscillation for a short-lived transient merged object. It is 
likely that peaks associated with QNMs appear after this peak, and the luminosity 
is eventually damped. For model (B-2), a long-lived transient object emits quasi- 
periodic gravitational waves for a certain duration. Consequently, the luminosity 
remains fairly large for a longer timescale than for model (C-2). 

In the simulations with models (B) and (C), the computation crashed soon after 
the formation of an apparent horizon, and hence it was not possible to observe 
gravitational waves associated with QNMs of the formed black hole. If it could be 
extracted, the peak associated with QNMs would appear at {t — rohs)/Pt=o ~ 2 for 
model (C-2) and 2.5 for model (B-2), before the luminosity is damped to zero. 

For the case of neutron star formation [model (A-2)], the evolution of the lumi- 
nosity is very different from that for the other models after the first peak is reached. 
In this case, the luminosity oscillates with average amplitude dE/dt ~ 5 x 10~^ after 
formation of the massive neutron star. We here estimate the approximate timescale 
for gravitational radiation reaction. 

The binding energy of the massive neutron star is approximately 

^~^ = 0.2m(^), ,4.7) 

where M is the gravitational mass and i? is a typical radius of the massive neutron 
star ~ 5M. The ratio of the kinetic energy T to is expected to be between 0.1 
and 0.2, so that 

T ~ (0.02 - 0.04)M 

Assuming that the energy luminosity remains ~ 10"*', as indicated in Fig. 10, and 
that the emission timescale is approximately derived by dividing T by dE/dt, we 
obtain 

Thus, within ^ 0.1 sec (or ^ 50P(=o), the massive neutron star will collapse into 
a black hole. As argued in Ref. 32), there are many other processes of dissipation 
and transport of angular momentum inside the massive neutron star, in addition to 
the emission of gravitational waves. However, the characteristic timescale for such 
processes is much longer than 1 sec. Thus, the emission timescale of gravitational 
waves would be shortest in this case. On the other hand, if tow were longer than 10 
seconds, other processes, such as magnetic braking effect, could be more important 
than the effect of gravitational radiation. 

All features mentioned above for F = 2.25 arc also found in gravitational wave- 
forms for r = 2. In Fig. 11, we display /i+ and hx for models (E-1) and (H-1). In 
Fig. 12, we also display the luminosity for these models. For model (H-1), a black 
hole is formed, while for (E-1), a neutron star is the outcome. As in the F = 2.25 
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case, this feature is well reflected in Figs. 11 and 12: Quasi-periodic oscillations are 
seen over a long duration for model (E-1), while for model (H-1), such oscillations 
are found only for a short timescale. 

For model (E-1), the luminosity remains ~ 10^^ after formation of a massive 
neutron star. Thus, the emission timescale of gravitational waves from a formed 
massive neutron star is likely shorter than 1 sec, as shown in Eq. (4-9). This implies 
that, as in the case of model (A), it eventually collapses into a black hole, due to the 
angular momentum dissipation through gravitational waves. 

4.2.2. Gravitational wave spectra 

In Fig. 13, spectra of gravitational waves are displayed for F = 2.25 [models 
(A-2), (B-2) and (C-2) on the left] and for T = 2 [models (F-1) and (H-1) on the 
right]. Since the computation crashed soon after formation of the black hole, the 
spectra of models (B-2), (C-2) and (H-1) are incomplete on the high frequency side 
[i.e., for /gw//qe ^ 4, where /qe is given in Eq. (4-1)]. In these cases, QNMs of 
black holes would be excited in the last stage of black hole formation. However, as 
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mentioned in §1, the frequency would be very high (/gw > 5 kHz ~ 5/qe)- Here, 
we focus only on the spectra for quasi-periodic oscillations of merged objects for 
which /gw ~ (2 — 3) kHz, so that the incompleteness with regard to the spectra 
on the higher frequency side is not a problem. For the spectra of model (A), it 
should be kept in mind that the peak associated with the quasi-periodic oscillations 
is underestimated in Fig. 13, because quasi-periodic waves will be emitted perhaps 
to t ~ 100Pf=o, as indicated in Eq. (4-9). For all models, we should also point out 



that in the real spectra of the coalescing binaries, /i+^x oc /gw^ /gw ^ /qe 
which is not taken into account in the present results. 

From Fig. 13, we find two typical frequencies of the quasi-periodic oscillation for 
each F: 
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/qpo ~ ' 



/qpo ~ ' 



1.8/qe « 1.7 kHz 



2.8/qe « 2.7 kHz 



1.8/,E^1.7kHzf2-«^«^ 
3.0/QE-2.9 kHzfl,'"^® 




for r = 2.25,(4-10) 



for r 



(4-11) 



This shows that the ratio of /qpo to /qe is not very sensitive to {M/R)oo. How- 
ever, this ratio for higher frequencies ~ 3/qe depends on the stiffness (or F) of the 
equations of state. The numerical results indicate that this ratio is larger for softer 
equations of state for a peak of higher frequency. (This fact agrees with the results 
of Newtonian and post Newtonian simulations. ^^-') As pointed out in Ref. 7), 
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from the frequencies of quasi-periodic oscillations, we could constrain the stiffness of 
equations of state. 

As explained above, the peaks of the Fourier spectra that arc associated with the 
quasi-periodic oscillations are determined by its accumulated cycles. Hence, if we 
would not find a large peak, we could conclude that a black hole is formed in a short 
timescale after merger. This implies that even without detecting QNMs or other 
signals of the black hole, it is possible to obtain evidence for black hole formation. 
On the other hand, if we would find a large peak, we could conclude that a massive 
neutron star is formed, at least temporarily. The criterion for the formation of 
black holes depends on the compactness and stiffness of the equations of state of 
the merging neutron stars. Thus, this information is also useful for constraining the 
equations of state. 

Prom gravitational waves emitted in the inspiraling stage with post-Newtonian 
templates of waveforms, the masses of the two neutron stars, and hence, the total 
mass, will be determined. ^'^^ As mentioned above, from the intensity of the peak of 
the spectra associated with quasi-periodic oscillations, we can infer the compactness 
(for a given stiffness of an equation of state). This combined observation could con- 
strain the relation between the compactness and mass of neutron stars. In this way, 
it is likely possible to strongly constrain nuclear equations of state using observed 
quantities such as the intensity of the Fourier peak of quasi-periodic oscillation, the 
ratio of /qpo to /qe, and the total mass of the system. 

Since /qpo is rather high, it may be difficult to detect quasi-periodic oscillation 
by first generation, kilometer-size laser interferometers, such as LIGO I. However, 
resonant-mass detectors and/or specially designed advanced interferometers, such as 
LIGO II, may be capable in the future to detect such high frequency gravitational 
waves. These detectors will provide us a wide variety of information on neutron star 
physics. 

4.2.3. Calibrations 

To assess the accuracy and robustness of our results as well as to investigate 
the influence of the approximate initial conditions on gravitational waveforms, we 
performed several test simulations. 

As mentioned above, the numerical accuracy of the gravitational wave amplitude 
depends on the location of the outer boundaries L and the grid spacing Ax. It is 
indicated that (1) smaller L{< A) and larger Ax result in the underestimation of /i+^x 
and (2) smaller L{<^ A) also induces spurious modulation of /i+,x (by which /i+^x 
deviate from zero systematically). However, we will show that with the typical grid 
size adopted in the present work (L ^ 0.6A), these numerical errors are sufficiently 
suppressed. 

First, we investigate the convergence of the wave amplitude resulting when the 
computational domain is enlarged. In Fig. 14 (a), we compare the merger waveforms 
for model (A) with (293,293,147) [(A-0)] and (505,505,253) [(A-2)] grid numbers. 
We note that L ^ Aqp for the smaller grid number and L ^ 2Aqp for the large 
grid number, where Aqp here denotes the wavelength of quasi-periodic gravitational 
waves emitted by a formed massive neutron star. Figure 14 (a) shows that the 
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Fig. 14. (a): /i+ and hy, for models (A-2) (solid curve) and (A-0) (dashed curve), (b): /i+ and hy, for 
models (C-3) (solid curve), (C-2) (long dashed curve) and (C-0) (dashed curve). The horizontal 
solid, dotted and dashed lines drawn near the vertical axis denote the theoretical amplitudes at 
t = computed in the 2.5, second post-Newtonian and quadrupole approximations. 



Table IV. Amplitudes of /i+,x at t = for several levels of the post-Newtonian approximation. 



Model 


quadrupole 


IPN 


1.5PN 


2PN 


2.5PN 


(A) 


0.755 


0.611 


0.738 


0.705 


0.681 


(B) 


0.766 


0.594 


0.760 


0.712 


0.676 


(C) 


0.785 


0.578 


0.793 


0.727 


0.671 


(E) 


0.735 


0.620 


0.712 


0.690 


0.675 


(H) 


0.759 


0.590 


0.752 


0.705 


0.670 



amplitude of gravitational waves from quasi-periodic oscillations is not sensitive to 
L. This indicates that a value of L that is ^ Aqp is sufficiently large for computing 
waveforms associated with quasi-periodic oscillations of a formed massive neutron 
star. On the other hand, it is not clear if convergence is achieved for the wave 
amplitude for t Pt=o- 

To evaluate the numerical error in the amplitude for t ^ Pt=o, we compare the 
waveforms for models (C-0), (C-2) and (C-3) in Fig. 14 (b), in which L/Xq = 0.36, 
0.63, and 0.79, respectively. In addition to the waveforms, the amplitudes predicted 
from the post-Newtonian theory are shown by horizontal lines. Here, the wave 
amplitude of the \m\ = 2 modes (sum of / = 2 to / = 6 modes) for a binary of equal 
point mass in circular orbits can be written as 



Mtotv 
D 



^ 17 2 o 3 15917 4 17 . 

1 + 2ttV^ TTV^ 

8 2880 4 



(4-12) 



at the 2.5 post-Newtonian order, ^^-^ *^ where Mtot is the sum of the gravitational 
masses of the two stars in isolation (hence Mtot 7^ -Madmo for a finite separation, due 

*-* In Ref. 35), post-Newtonian waveforms are given up to second order. To derive 2.5 post- 
Newtonian waveforms, we need to use the 2.5 post-Newtonian luminosity derived in Ref. 36). 
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Fig. 15. Left: h+ and as functions of the retarded time for models (A-2) (solid curve) and (A-3) 
(dashed curve). Right: The same as (a) but for (B-2) (solid curve) and (B-3) (dashed curve). 



to the existence of the binding energy between the two stars) and v = [Mtotfio)^/^ . 
We Ust the amphtudes of h+^x for several levels of the post-Newtonian approximation 
in Table IV. It is seen that the amplitude converges to ~ 0.7 for each model. Figure 14 
(b) indicates that for model (C-0), the amplitude is underestimated by 30-40 %, but 
as L increases, a convergence is achieved. The wave amplitudes for models (C-2) and 
(C-3) are nearly coincident and also agree with the 2.5 post-Newtonian amplitudes 
within < 10% error. Thus, with these computational domains with L ^ O.6A0, the 
amplitude seems to be computed within 10% error, and hence we may assume that 
the error on the wave amplitude calculated for models (B-2), (C-2), and (H-l), shown 
in Figs. 9 and 11, is less than 10%. This conclusion agrees with that of our recent 
numerical experiment on gravitational waveforms from inspiraling binary neutron 
stars in quasiequilibrium states. 

Next, we investigate the effect of the grid resolution {Ax) on gravitational wave- 
forms. In Fig. 15, we show /i+ and hx for models (A-2) and (A-3) and models (B-2) 
and (B-3) for comparison. At a glance, we find that the global properties in the wave- 
forms are similar with regard to results for high and low resolution. However, two 
differences are also found : (1) since the merger starts earlier, due to larger numeri- 
cal dissipation, merger waveforms start earlier for simulations with lower resolution; 
(2) the wave amplitude is underestimated for lower resolution. Property (2) can be 
understood in the following manner. In the simulation with lower resolution, the 
density gradient is less accurately computed. The wave amplitude becomes higher 
in association with the formation of higher density peaks. Thus, the lower resolution 
results in underestimation of the gravitational wave amplitude. [This property has 
been found in many Newtonian simulations using Eulerian codes; e.g., Ref. 9).] 

It is desirable to further investigate the convergence improving the resolution in 
order to draw definite conclusions about whether the convergence of the waveforms 
is fully achieved. However, it is pragmatically difficult to do this with the current 
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Fig. 16. Left: h+ and hx for models (A-0) (solid curve) and (A-1) (dashed curve). Right: The 
same as (a) but for models (B-0) (solid curve) and (B-1) (dashed curve). The amplitude is 
underestimated by 30-40 % in the early stages {t ^ Pt^o), because the computational domain 
in these simulations is too small (see Figs. 9 and 14 for comparison). 



computational resources, because a simulation with a large grid size, e.g., (633, 633, 
317), would take a very long computational time, perhaps ^ 200-250 CPU hours, 
which is about 1/3 of the computational time assigned to us. In the future, it will be 
feasible to improve the resolution using better computational resources. In addition, 
an adaptive mesh refinement (AMR) technique for high density regions might be 
necessary to save computational time and memory. ^'^^ 

Finally, we investigate the effect of an approaching velocity, which we ignored 
in our main simulations for F = 2.25. In this test, we perform a simulation with a 
large approaching velocity as fapp = O.lvorb for {M/R)oo = 0.12 [model (A-1)], and 
for iM/R)oo = 0.14 [model (B-1)] with (313,313,157) grid number. For model (B-1), 
we also decrease the reduction factor of the angular momentum by 1% (see Table 
II). In Fig. 16, we show and hx for these two cases. For comparison, we also plot 
the waveforms of models (A-0) and (B-0). Because of the approaching velocity, the 
waveforms before merger starts are very different for the two cases; the wavelength 
becomes short for models (A-1) and (B-1) soon after the simulations start. However, 
the quasi-periodic waveforms after the onset of the hydrodynamic interaction of 
two neutron stars are quite similar. This is because the quasi-periodic waves are 
associated with the fundamental oscillation modes of merged objects, which appear 
to depend weakly on the approaching velocity. Thus, we may conclude that the 
waveforms in the merger do not depend strongly on the approaching velocity, as 
long as Va_pp is not extremely large, i.e., as long as < O.lUorb- 

We note that /i+ in Fig. 16 deviates from zero in the late stages for t — Zobs ^ 
1.5Pt=o [in particular for model (B)]. This spurious effect results from the fact that 
the outer boundaries are too close to the strong field zone. Indeed, this numerical 
effect disappears when we increased the grid number to enlarge the computational 
domain (compare with Fig. 9). 
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§5. Summary 

Using a new large scale supercomputer at NAOJ, we have performed fully GR 
simulations for the merger of equal mass binary neutron stars, particularly focusing 
on gravitational waveforms. Since the computational domain is much larger than 
that used in previous simulations, gravitational waveforms can be computed much 
more accurately than those reported in Refs. 13) and 14). In the following, we 
summarize the results obtained in this paper. 

The merger process and final products depend strongly on the compactness and 
the stiffness of the equations of state for the merging neutron stars. For binaries of 
less compact neutron stars, massive neutron stars are formed, at least temporarily. 
In contrast, black holes arc formed in the merger of sufficiently compact neutron 
stars on a dynamical timcscale. However, the formation timescale of a black hole 
depends strongly on the compactness of the neutron stars. The criterion regarding 
compactness for prompt formation of a black hole depends also on the stiffness of 
equations of state. For F = 2.25 and 2, threshold values are ~ 0.16 and ~ 0.14, 
respectively. 

The nature of the gravitational waveforms during the merger depend sensitively 
on the compactness (or mass) of the merging neutron stars. For mergers of less 
compact binaries, such as in models (A), (E) and (F), a transient massive neutron 
star is formed, and survives for the emission timescale of gravitational radiation, 
which is much longer than the dynamical timescale. Such massive neutron stars are 
highly nonaxisymmctric, so that quasi-periodic gravitational waves associated with 
nonaxisymmetric deformation are emitted. These quasi-periodic waves have typical 
frequencies of ~ 2 — 3 kHz, which yield Fourier peaks in the frequency domain of 
gravitational waves. The ratio of the peak frequency to the frequency at innermost 
orbits, /qpo//qe5 is not sensitive to the compactness, but depends on the stiffness 
of the equations of state. This indicates that its observation could constrain the 
stiffness of the equations of state. It is also found that the luminosity of quasi- 
periodic gravitational waves is fairly large, so that a massive transient neutron star 
likely collapses into a black hole eventually, due to the angular momentum dissipation 
through gravitational radiation. 

In models of slightly larger compactness, such as models (B) and (G), a transient 
massive object is also formed after merger sets in. This massive object is also highly 
nonaxisymmetric. Since the lifetime of such a transient massive object is fairly long, 
it causes characteristic peaks associated with the quasi-periodic oscillation in the 
Fourier domain of gravitational waves. However, for models with large compactness, 
such as models (C), (D), (H) and (I), the merged object collapses into a black hole 
on a dynamical timescale (i.e., on a few rotational periods of the merged objects), 
and hence quasi-periodic gravitational waves are excited only on a short timescale. 

The intensity of the peaks in the Fourier spectra of gravitational waves associated 
with the nonaxisymmetric, quasi-periodic oscillations of the merged object depends 
on the lifetime of the transient massive object. If the transient object is long-lived, 
the peak becomes high, while if it collapses into a black hole in a few oscillation 
periods, the peaks are small. Therefore, from the intensity of the peak associated 
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with the quasi-periodic osciUation, it may be possible to determine the object formed 
after the merger in future observations. 

To this time, we have performed simulations focusing on binaries of two identical 
neutron stars. Although all the binary neutron stars observed so far are composed 
of neutron stars of two nearly identical masses, ^^-^ it seems that there is no fun- 
damental reason for the production of such symmetric systems. In the merger of 
two neutron stars of unequal mass, the merger process would be different from that 
for identical neutron stars. Since a less massive star is less compact, it is likely 
to be tidally disrupted before the dynamical instability of the orbital motion sets 
in. In this case, the tidally disrupted star could form accretion disks around the 
more massive companion. Sufficient accretion to the massive companion is likely 
to trigger gravitational collapse and the formation of a black hole. As a result, a 
system consisting of a black hole surrounded by accretion disks might be formed. 
Associated with the change of the merger process, gravitational waveforms would be 
modified.^") To survey possible outcomes in the merger of binary neutron stars, it is 
obviously necessary to perform simulations for two neutron stars of unequal mass. 

It is also desirable to improve the implementation of the initial conditions. In 
simulations carried out to this time, we have used quasiequilibrium states of a confor- 
mally flat three-metric as the initial conditions for simplicity. The conformal flatness 
approximation becomes a source of a certain systematic error when attempting to 
obtain realistic quasiequilibrium states. As a result, this approximation introduces 
a systematic error on the initial conditions and subsequent merger simulation. Since 
the magnitude of the ignored terms in the conformal flatness approximation seems 
to be small, we expect that this effect is not very serious. However, this conclusion 
is not entirely certain. To rule out this problem, it is necessary to perform sim- 
ulations using more realistic quasiequilibrium states of generic three geometries as 
initial conditions. 

In addition, there are technical issues requiring improvement. As discussed above 
and in the Appendix, black hole forming region does not have good resolution in our 
current computation. Obviously, it is necessary to improve the resolution around the 
black hole forming region. Since we have to prepare a large computational domain 
L, which is at least equal to the wavelength of gravitational waves, using restricted 
computational speed and memory, it is desirable to develop numerical techniques 
such as the AMR technique or nested grid technique to overcome this problem. This 
is also an issue to be resolved in the future. 
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Appendix A 

Accuracy 

We have monitored the degree of violation of the Hamiltonian constraint, conser- 
vation of the baryon rest-mass, ADM mass and angular momentum during numerical 
simulation. 

We have found that is less than 0.1 for a region in which p is larger than ~ 

Pmax) where pmax. IS the maximum value of p. For a less dense region, is often 
0(0. 1), because low density regions are not well resolved in the finite differencing 
scheme we used for the hydrodynamic equations. When the merged object starts 
collapsing, increases to 0(0. 1), even in a high density region. This indicates that 
the resolution in the black hole forming region is not very good either. 

The baryon rest-mass was found to be conserved within the truncation error 
in all computations. This is because no mass is ejected outside the computational 
domain, which is wide enough. 

In the case of massive neutron star formation, the ADM mass is conserved 
within 10% error throughout the simulation. At the last moment, at which excessive 
numerical error seems to accumulate, the simulation crashes, and the error of the 
ADM mass suddenly diverges, as shown in Fig. 17. Such a crash took place typically 
at some time in the range t = 700 - 900Madm (see §4). In the case of black hole 
formation, the ADM mass is conserved within 10% error until the merged object 
starts collapsing into a black hole. The error systematically increases whenever the 
lapse at the center approaches to zero, and the error in the conservation of the 
ADM mass is typically ~ 20% at the time of the formation of the apparent horizon 
(see Fig. 17). After this formation, the error rapidly increases to 0(1), and the 
computation crashes, at which time the error is of order unity. The main source 
of the error in the ADM mass appears to be the volume integral of Rj^^, because it 
includes the second spatial derivative of ^jj , and hence loses accuracy easily whenever 
7ij has a steep gradient. 

The angular momentum is conserved within 5% error throughout the simu- 
lation in the case of massive neutron star formation, as discussed in §4 (except for 
the last moment, at which excessive numerical error accumulates). For black hole 
formation cases, the angular momentum is conserved within ~ 5% error until the 
merged object starts collapsing. The error increases when the lapse at the center 
approaches to zero. The error for conservation in the angular momentum is typically 
~ 10% at the time of the formation of the apparent horizon. After the formation, 
the error increases to 0(1). 

The divergence of the numerical error after formation of the apparent horizon is 
likely due either to the lack of numerical resolution around black hole forming region 
or to the inappropriate choice of the spatial gauge. However, at present it is not 
clear to us whether one of these is the actual source. 
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